subroutine readinput
    use allocs
    use define
    USE IFPORT
    implicit none

    include "pspec.FI"
    include "files.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "accm.FI"
    include "mcbase.FI"
    include "FMC.FI"
    include "statistics.FI"
    include "GroupConstant.FI"
    include "mpi_arrays.FI"
!$  include "mtbuf.FI"
!$  include "omp_lib.h"
    integer :: unit, unitxsec, i, j, k, l, index, index_cmfd, g, gp, ig, ichar1st, temp, ratio, base, left, right
    integer :: ia, ja
!    character*200 logstring
    character*256 oneline, cardname
    character(1) probe, sc(1000)
    character(1) DOT,BANG,BLANK,SLASH,AST,POUND
    parameter (DOT='.',BANG='!',BLANK=' ',SLASH='/',AST='*')
    parameter (POUND='#')
    integer :: pin(100), id, val, npasm(100)
    integer :: nasm, ndataf
    integer :: totasm, totpin, rcp, asmrd, ntmpin, ntmasm
    real(8) :: pinsize, asmsize
    integer, parameter :: fref = 90
    integer :: i_c
    real(8) :: hr
    integer :: pxs, pxe, pys, pye, ix, iy
    integer :: yj, tcomp
    equivalence (probe,oneline)
    equivalence (sc(1),oneline)
    integer :: temprc(0:10)
    real(8) :: temprr(1:10)
    integer :: iis, iie
    integer :: hffx, hffy
    character*200 strseed
    integer :: stx, sty
    integer :: midpin
    real(8) :: nmod
   
    
    nskip=0;nprd=1
    
    localfn = ' '
    
!    call getarg(1, localfn)
!    ichar1st=ichar(localfn(1:1))
!    if(ichar1st.eq.0 .or. ichar1st.eq.32) localfn='input.txt'
    localfn='input.txt'
    
    unit =FIinp    
    unitxsec=unit+1
    open(unit, file=localfn,status='old', shared)

    read(unit,*) CMFDIACT, CMFDACT, nthrd

    read(unit,*) nax, nay, ndimpa
    read(unit,*) asmsize, nmpin, rcp
    read(unit,*) albedo(3:4), albedo(1:2)
    read(unit,*) history, inactive, active
    RSTDHAT = .false.
    ndelayFDB = 0
    fdbinact = CMFDIACT
    fdbact = CMFDACT
    refltype_asm = CENT

!$  call omp_set_dynamic(.FALSE.)
!$  call omp_set_num_threads(nthrd)
!$  write(logstring,  '(" Total # of Thread : ",i)') nthrd
    call log_(logstring)
    
    npx = ndimpa            ! # of pin div for a asm, x-dir
    npy = ndimpa            ! # of pin div for a asm, y-dir
    nda = ndimpa*ndimpa     ! # of total pin in a asm
    pinsize = asmsize/ndimpa
    hpx = pinsize           ! pin length
    hpy = pinsize           ! 
    vpin = hpx*hpy
    hax = pinsize * npx     ! asm length
    hay = pinsize * npy
    hcx = hax * nax         ! core length
    hcy = hay * nay
    hfmx = hpx/nmpin        ! fine mesh length        
    hfmy = hpy/nmpin     
    vfm  = hfmx*hfmy        ! vol for fine mesh
    nfmx = npx*nax*nmpin    ! # of fine mesh x-dir
    nfmy = npy*nay*nmpin
    rcx  = rcp              ! ratio fine to coarse x-dir
    rcy  = rcp
    hcmx = hfmx*rcx         ! coarse mesh length        
    hcmy = hfmy*rcy     
    vcm  = hcmx*hcmy        ! vol for coarse mesh
    ncmx = nfmx/rcx         ! # of coarse mesh x-dir 
    ncmy = nfmy/rcy
    totasm = nax*nay
    totpin = nax*npx*nay*npy
    a2m=ndimpa*nmpin
    
    npowpinref = nfmx*nfmy
    npowasmref = ncmx*ncmy

    allocate(core(nax, nay))
    allocate(assembly(npx, npy, 0:100))
    assembly(:,:,0)=0
    
   

    
    
    !if (SRCSPLT <1 .or. SRCSPLT>2) SRCSPLT = 1
    if (SRCSPLT <1) SRCSPLT = 1
    npasm = 0
    
    ! default CMFD convergence criteria
    max_out=1000; max_in=20;
    err_out=1e-5_8; err_in=1e-2_8;
    ntypepin = 0
    nmxr = 0
    ireft = 0
    WFSDinact=0; WFSDact=0; Woffset=0
    lnoacc = .false.
    rungcmfd = .false.
    hybridcmfd = .false.
	lreadFDM = .false.
	nignracc = 0
	mltfinact = 1.0
	ppreconst = .false.
	CMFDDISCARD = .false.
	WBINCMFD = .false.
	WERRDIST = .false.
	
	ntypeasm = 0
	
	MC_ASMCAL=.false.
	
    allocate(typepin(MAX_TYPEPIN))
    do while (TRUE)
        read(unit,'(a256)') oneline
        if(probe.eq.DOT .or. probe.eq.SLASH) exit
        if(probe.eq.BANG .or. oneline.eq.BLANK .or. probe.eq.POUND) cycle
        
        read(oneline, *) cardname
        call toupper1(cardname)
        select case(cardname)
            case('CMFDCONV')
                read(oneline, *) cardname, max_out, err_out, max_in, err_in
            case('FDBCOND')
                read(oneline, *) cardname, nskip, nprd
            case('PIN')
                read(oneline, *) cardname, id, typepin(id)%nr, temprr(1:typepin(id)%nr), temprc(0:typepin(id)%nr)
                forall(i=0:typepin(id)%nr) typepin(id)%rc(i)=temprc(typepin(id)%nr-i)
                forall(i=1:typepin(id)%nr) typepin(id)%rr(i)=temprr(typepin(id)%nr-i+1)
                if (nmxr<typepin(id)%nr) nmxr=typepin(id)%nr
                do i=typepin(id)%nr,1,-1
                    typepin(id)%rv(i)=PI*typepin(id)%rr(i)**2
                end do
                typepin(id)%rv(i)=hpx*hpy
                do i=0,typepin(id)%nr-1
                    typepin(id)%rv(i)=typepin(id)%rv(i)-typepin(id)%rv(i+1)
                    typepin(id)%rvr(i)=1._8/typepin(id)%rv(i)
                enddo
                typepin(id)%rvr(i)=1._8/typepin(id)%rv(i)
                ntypepin=ntypepin+1
                !pin(id) = val
            case('ASSEMBLY')
                id = 0
                read(oneline, *) cardname, id
                if (id .eq. 0) then
                    write(*,*) "Error in Assembly Index!"
                    exit 
                endif
                do j=1, npy
                    read(unit,'(a256)') oneline
                    if(probe.eq.DOT .or. probe.eq.SLASH) exit
                    if(probe.eq.BANG .or. oneline.eq.BLANK .or. probe.eq.POUND) cycle
                    ntmpin=nfields(oneline)
                    if (j .eq. 1 .and. ntmpin .eq. nda) then
                        read(oneline,*) assembly(:,:,id)
                        exit
                    elseif (ntmpin .ne. npx) then
                        write(*,*) "Error in # of pins in an assembly!", id
                        exit
                    endif
                    read(oneline, *) assembly(:,npy-j+1,id)
                enddo
                if (id>ntypeasm) ntypeasm =id
!                if (npasm .ne. 0 .and. npasm .ne. ndataf) write(*,*) "Geometry error : assembly ", id
!                npasm = ndataf
            case('CORE')
                asmrd = 0
                do j=1, nay
                    read(unit,'(a256)') oneline
                    if(probe.eq.DOT .or. probe.eq.SLASH) exit
                    if(probe.eq.BANG .or. oneline.eq.BLANK .or. probe.eq.POUND) cycle
                    ntmasm=nfields(oneline)
                    if (j .eq. 1 .and. ntmasm .eq. totasm) then
                        read(oneline,*) core(:,:)
                        asmrd = ntmasm
                        exit
                    elseif (ntmasm .ne. nax) then
                        write(*,*) "Error in # of assemblies in the core!"
                        exit
                    endif
                    read(oneline, *) core(:,nay-j+1)
                    asmrd = asmrd  + nax
                enddo
                if (totasm .ne. asmrd) then
                    write(*,*) "Not enough assemblies in the core."
                    cycle
                endif
            case("MCASMCAL")
                MC_ASMCAL=.true.
                homtype=HT_PTASM
                dopprec_fdm=.false.
                do while (TRUE)
                    read(unit,'(a256)') oneline
                    if(probe.eq.DOT .or. probe.eq.SLASH) exit
                    if(probe.eq.BANG .or. oneline.eq.BLANK .or. probe.eq.POUND) cycle
                    
                    read(oneline, *) cardname
                    ntmpin=nfields(oneline)
                    call toupper1(cardname)
                    select case(cardname)
                    case("HOMTYPE") 
                        read(oneline, *) cardname, homtype
                    case("GROUPS") 
                        allocate(gtog(ntmpin-1))
                        read(oneline, *) cardname, gtog(1:ntmpin-1)
                    case("MCOPT")
                        read(oneline, *) cardname, nht_asm, inactcyc_asm, actcyc_asm
                    case("FDM_PPREC")
                        read(oneline, *) cardname, dopprec_fdm
                    endselect
                enddo
                if(homtype .eq. HT_PTASM)  dopprec_fdm=.false.

            case ("CMFDWLDT")
                read(oneline, *) cardname, lwldt
            case ("RESETCMFD")
                RSTDHAT = .true.
            case ("DELAYEDFDB")
                read(oneline, *) cardname, ndelayFDB
            case ("CMFDFEEDBACK")
                read(oneline, *) cardname, fdbinact, fdbact
            case ("UNIFORMFSRC")
                print *, "[ Forced uniform fission source neutron ]"
                unifsrc = 1
            case ("RESAMPLING")
                lresampling = .true.
            case ("MGFDMFIRST")
                lmgfdmfirst = .true.
            case ("TALLYDISCARD")
                print *, "[ Tallies are discarded from accumulation for CMFD ]"
                lnoacc = .true.
            case ("FSDCYCLE") 
                read(oneline, *) cardname, WFSDinact, WFSDact, Woffset, WFSDacc
            case ("NIGNRACC") 
                read(oneline, *) cardname, nignracc
            case ("RUNGCMFD")
                rungcmfd = .true.
                print *, "[ GCMFD is performed instead of CMFD ]"
            case ("HYBRID")
                hybridcmfd = .true.
			case ("READFDM")
				lreadFDM = .true.
			case ("REFLECTEDGE")
			    refltype_asm = EDGE
			case ("MLTFINACT")
			    read(oneline, *) cardname, mltfinact
			case ("PPRECONST")
			    read(oneline, *) cardname, nskip_p2b
			    ppreconst = .true.
			case ("CMFDDISCARD")
			    CMFDDISCARD=.true.
			    read(oneline, *) cardname, ndscrd
			case ("WBINCMFD")
			    WBINCMFD=.true.
			case ("WERRDIST")
			    WERRDIST=.true.
            case ("REFERENCE")
                read(oneline, *) cardname, ireft
                if (ireft .eq. 1) then
                    call dmalloc(refp2bf, nfmx, nfmy)
                    open(2001, file='REF_P2BF.bin', form='binary', status='old')
                    read(2001) refp2bf
                    close(2001)
                    sump2bf=sum(refp2bf)
                    call dmalloc(refpsi, nfmx, nfmy)
                    open(2001, file='REF.bin', form='binary', status='old')
                    read(2001) refpsi
                    close(2001)
                    pinpowsum = sum(refpsi)
                    refpsi = refpsi/pinpowsum
                    npowpinref = 0
                    do j=1, nfmy
                        do i=1, nfmx
                            if(refpsi(i,j)>0) npowpinref = npowpinref+1
                        enddo
                    enddo

                    call dmalloc(refcmfd, ncmx*ncmy)
                    l=0
                    do j=1, ncmy
                        do i=1, ncmx
                            l=l+1
                            stx=(i-1)*rcx
                            sty=(j-1)*rcy
                            refcmfd(l)=sum(refpsi(stx+1:stx+rcx, sty+1:sty+rcy))
                        enddo
                    enddo
                    refcmfd=refcmfd/sum(refcmfd)
                    call dmalloc(powasm, nax, nay)
                    
!                    open(2001, file='REF_CMFD.bin', form='binary', status='old')
!                    read(2001) trefcmfd
!                    close(2001)
!                    pinpowsum = sum(trefcmfd)
!                    trefcmfd = trefcmfd/pinpowsum
!                    npowasmref = 0
!                    do j=1, ncmy
!                        do i=1, ncmx
!                            if(trefcmfd(i, j)>0) npowasmref = npowasmref+1
!                        enddo
!                    enddo
                    
                endif
        end select
    enddo
    close(unit)
    
    ! treatment for the reflection when the axis is on the assemlby center============================================================        
    call dmalloc(volfr, nfmx, nfmy)
    allocate (volfr_r( nfmx, nfmy))
    volfr=1
    volfr_r=1._8
    ibndry(DW)=1
    ibndry(DE)=nfmx
    ibndry(DN)=nfmy
    ibndry(DS)=1
    if (refltype_asm .eq. CENT)then
        if (albedo(DW) .eq. 0._8) ibndry(DW) = ndimpa/2+1
        if (albedo(DN) .eq. 0._8) ibndry(DN) = nfmy-ndimpa/2-1
        refltype_pin=EDGE
        if (mod(ndimpa,2) .eq. 1) then
            refltype_pin=CENT
            midpin = ndimpa/2+1
            allocate(reflpin(nfmx, nfmy))
            reflpin =0
            if (albedo(DW) .eq. 0._8) then
                reflpin(midpin,:) =1
                volfr(midpin,:) =volfr(midpin,:)*0.5
                volfr_r(midpin,:) =volfr_r(midpin,:)*2
            endif
            midpin = nfmy-ndimpa/2
            if (albedo(DN) .eq. 0._8) then
                reflpin(:,midpin) =reflpin(:,midpin)+2
                volfr(:,midpin) =volfr(:,midpin)*0.5
                volfr_r(:,midpin) =volfr_r(:,midpin)*2
            endif
        endif
    endif
    !=================================================================================================================================        
    write(logstring, '(" # of fine mesh    : ",i)') nfmx*nfmy
    call log_(logstring)
    write(logstring, '(" # of coarse mesh  : ",i)') ncmx*ncmy
    call log_(logstring)
    
    sizeprb = vcm*ncmx*ncmy
    write(logstring, '(" size of problem   : ",f12.2)') sizeprb 
    call log_(logstring)

!    do i=1, totasmdksl
!        dh_cmfd(i) = pinsize*  npasm(core(i))
!    enddo
    
    ! prepare map for tilt check ===================
    hffx=nfmx/2; hffy=nfmy/2
    call dmalloc(qtmapx, nfmx)
    call dmalloc(qtmapy, nfmy)
    qtmapx(1:hffx)=1
    qtmapx(nfmx-hffx+1:nfmx)=2
    qtmapy(1:hffy)=1
    qtmapy(nfmy-hffy+1:nfmy)=2
    call dmalloc0(fis1_4,0,2,0,2,1,nthrd)
    call dmalloc0(fis1_4acc,0,2,0,2,1,nthrd)
    !===============================================
    
    call dmalloc(pin2bf, nfmx, nfmy)
    call dmalloc (cplst,   nfmx, nfmy)
    call dmalloc0(hacx, 0, nfmx)
    call dmalloc0(hacy, 0, nfmy)
    

!    call dmalloc (compl_cmfd, totm_cmfd)
!    call dmalloc (hcmx, ncmx)
!    call dmalloc (hcmy, ncmy)
!    call dmalloc0(hac_cmfdx, 0, ncmx)
!    call dmalloc0(hac_cmfdy, 0, ncmy)
    
    call dmalloc(psif,   nfmx, nfmy)
    call dmalloc(psifp,  nfmx, nfmy)
!    call dmalloc(psia, nfmx*nfmy,1)

    call dmalloc(ppf,    nfmx, nfmy)
    call dmalloc(srcdst, nfmx, nfmy)
    call dmalloc(adjp,   ncmx, ncmy)
	call dmalloc(initfsd, nfmx, nfmy)
    adjp=1._8
    call dmalloc(fiswtmap, ncmx, ncmy)
    fiswtmap=1    
   
    index = 1
    hacx(0) = 0.
    hacy(0) = 0.
    
    do j=1, nay
        do i=1, nax
            do iy=1, npy
                pys = (j-1)*npy*nmpin + nmpin*(iy-1)+1
                pye = (j-1)*npy*nmpin + nmpin*iy
                do ix =1, npx
                    pxs = (i-1)*npx*nmpin + nmpin*(ix-1)+1
                    pxe = (i-1)*npx*nmpin + nmpin*ix
                    cplst(pxs:pxe, pys:pye) = assembly(ix, iy, core(i, j)) 
                end do
            end do
        end do 
    end do
    
    
    hr = 1._8/(pinsize/nmpin)
    hacx = 0
    do i=1, nfmx
        hacx(i) = hacx(i-1)+hfmx
    end do
    hacy = 0
    do i=1, nfmy
        hacy(i) = hacy(i-1)+hfmy
    end do    
    

    
    unitxsec =FXsec
    open(unitxsec, file='xsec.txt',status='old')

    read(unitxsec,'(a256)') oneline
    if (nfields(oneline) .eq. 3) then
        read(oneline,*) ng, nComp, mge
    else
        read(oneline,*) ng, nComp
        mge=ng;
    endif
    
    
    ! MGCMFD test %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    lmgcmfd = .false.
    ngcmfd = 2
    ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    call dmalloc(gc2, ng)
    call dmalloc(gcmc, ng)
    gc2(1:mge)=1
    gc2(mge+1:ng)=2
    if (ngcmfd .eq. 2)then
        gcmc=gc2
    elseif (ngcmfd .eq. ng) then
        do i=1, ng
            gcmc(i)=i
        enddo
    else
    endif
!    else
!        do i=1, ng
!            gc2(i)=i
!        enddo
!    endif
    
    call dmalloc(autocor_arr, 3, inactive+active)

#ifdef _THFDB
    call dmalloc(Xdif, totpin, 2)
    call dmalloc(Xrmv, totpin, 2)
    call dmalloc(Xnuf, totpin, 2)
    call dmalloc(Xchi, totpin, 2)
    call dmalloc(Xsct, totpin, 2)
    call dmalloc(Xtot, totpin, 2)
    call dmalloc(Xabs, totpin, 2)
    call dmalloc(XsctM, 2, 2, totpin)
    call dmalloc(cplstth,nfmx,nfmy)
    ! assign assembly index to composition
!    do j=1, nfmy
!        ja = int4((j-1)/a2m)+1
!        do i=1, nfmx
!            ia=int4((i-1)/a2m)+1
!            cplst(i,j)=(ja-1)*nay+ia
!        enddo
!    enddo
    tcomp=1
    do j=1, nfmy
        do i=1, nfmx
            cplstth(i,j)=cplst(i,j)
            cplst(i,j)=tcomp
            tcomp=tcomp+1
        enddo
    enddo

#else
    call dmalloc(Xdif,  nComp, ng)
    call dmalloc(Xrmv,  nComp, ng)
    call dmalloc(Xnuf,  nComp, ng)
    call dmalloc(Xchi,  nComp, ng)
    call dmalloc(Xsct,  nComp, ng)
    call dmalloc(Xtot,  nComp, ng)
    call dmalloc(Xtotr, nComp, ng)
    call dmalloc(Xnfot, nComp, ng)
    call dmalloc(Xabs,  nComp, ng)
    call dmalloc(Xkappa,nComp, ng)
    call dmalloc(XsctM, ng, ng, nComp)
    do i=1 , nComp
        read(unitxsec,*)
        do g=1, ng
            read(unitxsec,*) temp, Xdif(i,g), Xrmv(i,g), Xnuf(i,g), Xkappa(i,g), Xchi(i,g)
        end do
        do g=1, ng
            read(unitxsec,*) temp, XsctM(1:ng, g, i)
            do gp=1, ng
                if (XsctM(gp, g, i) < 0._8) XsctM(gp, g, i)=0._8
            enddo
            XsctM(g,g,i) = 1._8/3/Xdif(i,g)-Xrmv(i,g)
            Xsct (i, g) = sum(XsctM(1:ng, g, i))
        end do
    end do
    close(unitxsec)
    Xtot=1._8/(3*Xdif)
    Xtotr=1._8/Xtot
    Xnfot=Xnuf/Xtot
    Xabs=Xtot-Xsct
    call dmalloc(cdfsct,ng,ng,nComp)
    do i=1, nComp
        do g=1, ng
            cdfsct(1, g, i) = XsctM(1, g, i)
            do ig=2, ng
                cdfsct(ig, g, i) = cdfsct(ig-1, g, i)+XsctM(ig, g, i)
            enddo
            cdfsct(:,g,i)=cdfsct(:,g,i)/cdfsct(ng,g,i)
        enddo
    enddo
    
#endif
    
    call dmalloc0(jf,       1,4, 0,nfmx+1, 0,nfmy+1, 1,ng)
    call dmalloc0(phif,     0,MAX_RING,1,nfmx,1,nfmy,1,ng)
    call dmalloc(phi_fnm,  nfmx,nfmy,  ng)
    
    call dmalloc(srcPDF,   ncmx,ncmy)

    call dmalloc0(trk_est,  0,MAX_RING,1,nfmx,1,nfmy,1,ng)
    

!! Variables for DF
!    call dmalloc(coef_df, totm_cmfd*2, ng)

    cntacc = 0
#ifdef MPIF 
    nmod = 1.0
    if (mltfinact> 1.0) nmod = mltfinact   
    
    allocate(ntrn(history*mltfinact*5/nprocs));
    allocate(nbuf(history*mltfinact/nprocs))
	if(myrank .eq. 0) print *, "NTRN", (sizeof(ntrn)+sizeof(nbuf))*1e-6_8, "mb"
#else
    allocate(ntrn(history*mltfinact*5));
	print *, "NTRN", sizeof(ntrn)*1e-6_8, "mb"

#endif
    phif = 1.    
    
   
    call dmalloc0(f2cx, 0,nfmx+1)
    call dmalloc0(f2cy, 0,nfmy+1)
    call dmalloc0(c2fx, 0,ncmx+1)
    call dmalloc0(c2fy, 0,ncmy+1)

    index =1
    index_cmfd = 1
    do i=1, ncmx
        c2fx(i) = index
        f2cx(index: index+ rcx-1) = i
        index = index+ rcx
    end do
    index =1
    index_cmfd = 1
    do i=1, ncmy
        c2fy(i) = index
        f2cy(index: index+ rcy-1) = i
        index = index+ rcy
    end do


    ! map for vacuum geometry 2010.03.30 
    call dmalloc(map_exp,ncmx*ncmy)
    call dmalloc0(map_exp2,0,ncmx+1,0,ncmy+1)
    ! -1 for vacuum assembly
    map_exp2(1:ncmx, 1:ncmy) = -1 
    call dmalloc(map_srk,ncmx*ncmy)
    call dmalloc(nplane, ncmy)
    iie=1;iis=1;
    do j=1, ncmy
        do i=1, ncmx
            if (cplst(c2fx(i), c2fy(j))>0) then
                map_exp(iie)=iis
                map_exp2(i,j)=iis
                map_srk(iis)=iie
                nplane(j)=nplane(j)+1 
!                if (ireft .eq. 1) refcmfd(iis)=trefcmfd(i,j)
                iis=iis+1
            endif
            iie=iie+1
        enddo
        nreal = nreal + nplane(j)
    enddo

!   use trefcmfd for assembly power check

!    if (ireft .eq. 1) deallocate(trefcmfd)
    
    
!    call dmalloc(jcvt2cmfdx, nfmx*2)
!    call dmalloc(jcvt2cmfdy, nfmy*2)
    call dmalloc(jf2jc, 4, nfmx, nfmy)
!    allocate(jf2jc(4, nfmx, nfmy))
    jf2jc = FALSE
    jf2jc(DN, :, nfmy) = TRUE
    jf2jc(DS, :, 1) = TRUE
    jf2jc(DW, 1, :) = TRUE
    jf2jc(DE, nfmx, :) = TRUE
    do j=1, nfmy
        do i=1, nfmx
            if (i .ne. nfmx )then
                if(f2cx(i) .ne. f2cx(i+1)) jf2jc(DE, i, j) = TRUE
            endif
            if (i .ne. 1 ) then
                if(f2cx(i-1) .ne. f2cx(i)) jf2jc(DW, i, j) = TRUE
            endif
            if (j .ne. nfmy ) then
                if(f2cy(j) .ne. f2cy(j+1)) jf2jc(DN, i, j) = TRUE
            endif
            if (j .ne. 1 ) then
                if(f2cy(j-1) .ne. f2cy(j)) jf2jc(DS, i, j) = TRUE
            endif
        end do
    end do
    
    
    
!    jcvt2cmfdx = FALSE
!    jcvt2cmfdy = FALSE
!    left = 1
!    do i=1, nfmx
!        right = left + rcx*2-1
!        jcvt2cmfdx(left ) = TRUE
!        jcvt2cmfdx(right) = TRUE
!        left = right+1
!    end do
!    left = 1
!    do i=1, nfmy
!        right = left + rcy*2-1
!        jcvt2cmfdy(left ) = TRUE
!        jcvt2cmfdy(right) = TRUE
!        left = right+1
!    end do
 
!    ! for test MC
!    call dmalloc(absn, totmesh, ng)
!    call dmalloc(srcn, totmesh, ng)
!    call dmalloc(sctn, ng, ng, totmesh)
    


!    if (PDFREF) then
!        logstring = "Read source distribution from reference solution"
!        call log_(logstring)    
!        open(fref, file='REF.txt', status='old')
!        call dmalloc(refphi, totmesh, ng)
!        do g = 1, ng
!            read(fref, *) refphi(:,g)
!        end do
!        close(fref)
!
!        call dmalloc(refc, totm_cmfd, ng)
!        refc = 0
!        do g=1, ng
!            do i=1, totmesh
!                i_c = fi2co(i)
!                ! cycle flux are used to obtain initial guess
!                refc(i_c, g) = refc(i_c, g) + refphi(i,g)*h(i)
!            end do
!            do i=1, totm_cmfd
!                refc(i, g) = refc(i, g)/h_cmfd(i)
!            end do
!        end do
!    endif   
    
!#ifdef _REFPSI
!    open(1985, file='refpsi.txt',status='old')
!    read(1985,*) ndrefpsi
!    call dmalloc(refpsi,ndrefpsi,ndrefpsi)
!    l2refpsi = 0._8
!    normrefpsi = 0._8
!    do j=1, ndrefpsi
!        yj = ndrefpsi+1-j
!        read(1985, *) refpsi(1:ndrefpsi, yj)
!        do i=1, ndrefpsi
!            l2refpsi = l2refpsi + refpsi(i,yj)**2
!            normrefpsi = normrefpsi + refpsi(i,yj)
!        enddo
!    enddo
!    l2refpsi = sqrt(l2refpsi)
!    
!#endif
    ! homogenization
    call dmalloc(hpsi ,ncmx,ncmy)
    call dmalloc(hphi ,2,ncmx,ncmy)
    call dmalloc(hj   ,4,2,ncmx,ncmy)
    call dmalloc(hDif ,2,ncmx,ncmy)
    call dmalloc(hTot ,2,ncmx,ncmy)
    call dmalloc(hRmv ,2,ncmx,ncmy)
    call dmalloc(hNuF ,2,ncmx,ncmy)
    call dmalloc(hmSct,2,2,ncmx,ncmy)
    call dmalloc0(hjac, 1,4, 1,2, 0,ncmx+1, 0,ncmy+1)


    call dmalloc(rphids ,2,ncmx,ncmy)
    call dmalloc(rtotds ,2,ncmx,ncmy)
    call dmalloc(rNuFds ,2,ncmx,ncmy)
    call dmalloc(rmsctds,2,2,ncmx,ncmy)
    call dmalloc0(rjds, 1,4, 1,2, 0,ncmx+1, 0,ncmy+1)

!    ! homogenized reaction for omp
!    call dmalloc(hphio,2,ncmx,ncmy,nthrd)
!    call dmalloc(hToto,2,ncmx,ncmy,nthrd)
!    call dmalloc(hNufo,2,ncmx,ncmy,nthrd)
!    call dmalloc(hmScto,2,2,ncmx,ncmy,nthrd)

    
    ! accumulation
    call dmalloc(psifAC,nfmx,nfmy)
!    call dmalloc(phivac,nfmx,nfmy,ng)
    call dmalloc0(phivac,0,nmxr,1,ng,1,nfmx,1,nfmy)
    call dmalloc(jfAC,4,nfmx,nfmy,ng)
    phivac = 1.
    
    ! statistics
    call dmalloc(stphi,nfmx,nfmy,ng)
    call dmalloc(stphisq,nfmx,nfmy,ng)
    call dmalloc(stpsiz, nfmy)
    call dmalloc(stpsizsq, nfmy)
    call dmalloc(stpsi,nfmx,nfmy)
    call dmalloc(stpsisq,nfmx,nfmy)
    call dmalloc(stpsia, nfmx*nfmy, 1)
    call dmalloc(stpsiasq,nfmx*nfmy,1)
    ! pin power
    call dmalloc(stpp,nfmx,nfmy)
    call dmalloc(stppsq,nfmx,nfmy)
    
    call dmalloc(trkdense, nfmy)
    
    
    call dmalloc(errdist, npowpinref)
    
!!    allocate(phif2AC(nfmx, nfmy, ng))
!    allocate(phisaAC(nfmx+1, nfmy+1, ng))
!    allocate(mmt1stAC(nfmx+1, nfmy+1, ng))
!    allocate(mmt2ndAC(nfmx+1, nfmy+1, ng))
!    allocate(mmtAbsrAC(nfmx+1, nfmy+1, ng))
!    allocate(mmtNuFAC(nfmx+1, nfmy+1, ng))
    call dmalloc(pxx,nfmx+1,nfmy+1)
    call dmalloc(pyy,nfmx+1,nfmy+1)
    call dmalloc(pxy,nfmx+1,nfmy+1, 2)
    call dmalloc(tlfmc, 12, 4, nfmx, nfmy)
    call dmalloc(ffmc, 12, 4, nfmx, nfmy) 
    call dmalloc(phi0, nfmx+1, nfmy+1)
!    call dmalloc(tpxx, nfmx+1, nfmy+1)
!    call dmalloc(tpyy, nfmx+1, nfmy+1)
    call dmalloc(tphi2, 2, nfmx+1, nfmy+1)
    call dmalloc(tj, 2, nfmx+1, nfmy+1)
    call dmalloc(txsurf, 2, nfmx+1, nfmy+1)
    call dmalloc(tAbs, nfmx+1, nfmy+1)
    call dmalloc(tNuf, nfmx+1, nfmy+1)
    call dmalloc(txnavg, nfmx+1, nfmy+1)
!    deallocate(assembly)

!$  call dmalloc0(trkomp ,0,nmxr ,1,ng ,1,nfmx, 1,nfmy, 1,nthrd)
!$  call dmalloc(psiomp  ,nfmx, nfmy, nthrd)
!$  call dmalloc(jfomp,4,ng,nfmx,nfmy,nthrd)
!$  call dmalloc(gkomp, nthrd)
#ifdef MPIF
!$  allocate(ntrnomp(history*5/nthrd/nprocs,2,nthrd))
#else
!$  allocate(ntrnomp(history*5/nthrd,2,nthrd))
#endif
!$  call dmalloc(nqomp,3,nthrd)

    
#ifdef MPIF
    ! MPI BUFFER
    call dmalloc(psifrcv, nfmx, nfmy)
    call dmalloc(totrcv, ngcmfd, ncmx, ncmy)
    call dmalloc(nufrcv, ngcmfd, ncmx, ncmy)
    call dmalloc(phivrcv, ngcmfd, ncmx, ncmy)
    call dmalloc(sctrcv, ngcmfd, ngcmfd, ncmx, ncmy)
    !call dmalloc(jcrcv, 4, ngcmfd, ncmx, ncmy)
	call dmalloc0(jcrcv,1,4, 1,2, 0,ncmx+1, 0,ncmy+1)
#else
#endif
    
    ! Reaction Rate Tally
!    call dmalloc(rrtot,  2,   ncmx, ncmy)
!    call dmalloc(rrnuf,  2,   ncmx, ncmy)
!    call dmalloc(rrabs,  2,   ncmx, ncmy)
!    call dmalloc(rrsct,  2,2, ncmx, ncmy)
!    call dmalloc(hmphi,  2,   ncmx, ncmy)
!    call dmalloc(hmphit, 2,   ncmx, ncmy)
!    call dmalloc(hmjc,   4,2, ncmx, ncmy)
!    call dmalloc(hmjcac, 4,2, ncmx, ncmy)
    call dmalloc0(jcomp, 1,4, 1,ngcmfd, 0,ncmx+1, 0,ncmy+1, 1,nthrd)
    
    
    ! prepare arrays for generating group constants
    
    
    ! allocate group constants
    if (lmgcmfd) then
        call dmalloc(GCTphiv,ng,ncmx,ncmy,nthrd)
        call dmalloc(GCTtot,ng,ncmx,ncmy,nthrd)
        call dmalloc(GCTnuf,ng,ncmx,ncmy,nthrd)
        call dmalloc(GCTabsr,ng,ncmx,ncmy,nthrd)
        call dmalloc(GCTsct,ng,ng,ncmx,ncmy,nthrd)
    else
        call dmalloc(GCTphiv,  2,ncmx,ncmy,nthrd)
        call dmalloc(GCTtot,   2,ncmx,ncmy,nthrd)
        call dmalloc(GCTnuf,   2,ncmx,ncmy,nthrd)
        call dmalloc(GCTabsr,  2,ncmx,ncmy,nthrd)
        call dmalloc(GCTsct, 2,2,ncmx,ncmy,nthrd)
    endif
    
    call dmalloc(GCAphiv,  2,ncmx,ncmy)
    call dmalloc(GCAtot,   2,ncmx,ncmy)
    call dmalloc(GCAnuf,   2,ncmx,ncmy)
    call dmalloc(GCAabsr,  2,ncmx,ncmy)
    call dmalloc(GCAsct, 2,2,ncmx,ncmy)
    
    ! for 2g cmfd
    call dmalloc(hj2,    4,2,ncmx,ncmy)
    call dmalloc(hphi2,    2,ncmx,ncmy)
    call dmalloc(htot2,    2,ncmx,ncmy)
    call dmalloc(hdif2,    2,ncmx,ncmy)
    call dmalloc(hrmv2,    2,ncmx,ncmy)
    call dmalloc(hnuf2,    2,ncmx,ncmy)
    call dmalloc(hmsct2, 2,2,ncmx,ncmy)

    ! mg->2g scattering matrix
    call dmalloc(GCsctm2g, 2,ng, ncomp)
    do i=1, ncomp
        do g=1, ng
            do gp=1, ng
                GCsctm2g(gc2(gp), g, i)=GCsctm2g(gc2(gp),g,i)+XsctM(gp,g,i);
            enddo
        enddo
    enddo
    hmzt = HMZMC
    
    ! arrays for Larsen's CMFD
    call dmalloc (phiau, 2, ncmy, ncmx)
    call dmalloc0(phisuu, 1,2, 0,ncmy, 0,ncmx)
    call dmalloc (tphiau, 2, ncmy, ncmx, nthrd)
    call dmalloc (tphisuu, 2, ncmy+1, ncmx+1, nthrd)
    
    ! arrays for GCMFD
    call dmalloc (trcwx,     2, ncmx,     ncmy)
    call dmalloc (trcwy,     2, ncmx,     ncmy)
    call dmalloc0(phi2x,   1,2, 0,ncmx, 0,ncmy)
    call dmalloc0(phi2y,   1,2, 0,ncmx, 0,ncmy)
    call dmalloc0(phi2trx, 1,2, 0,ncmx, 0,ncmy)
    call dmalloc0(phi2try, 1,2, 0,ncmx, 0,ncmy)
    ! tallies for GCMFD
    call dmalloc (ttrcwx,    2, ncmx,     ncmy,   nthrd)
    call dmalloc (ttrcwy,    2, ncmx,     ncmy,   nthrd)
    call dmalloc0(tphi2x,  1,2, 0,ncmx, 0,ncmy, 1,nthrd)
    call dmalloc0(tphi2y,  1,2, 0,ncmx, 0,ncmy, 1,nthrd)
    call dmalloc0(tphi2trx,1,2, 0,ncmx, 0,ncmy, 1,nthrd)
    call dmalloc0(tphi2try,1,2, 0,ncmx, 0,ncmy, 1,nthrd)
    ! MPI arrays for  GCMFD
    call dmalloc (trcwxrcv,    2, ncmx,     ncmy)
    call dmalloc (trcwyrcv,    2, ncmx,     ncmy)
    call dmalloc0(phi2xrcv,  1,2, 0,ncmx, 0,ncmy)
    call dmalloc0(phi2yrcv,  1,2, 0,ncmx, 0,ncmy)
    call dmalloc0(phi2trxrcv,1,2, 0,ncmx, 0,ncmy)
    call dmalloc0(phi2tryrcv,1,2, 0,ncmx, 0,ncmy)
    
    ! tallies for assembly decoupled parallelization
#ifdef ASMPRL    
    call dmalloc(tnfs, ncmx, ncmy, nthrd)
    call dmalloc(tnin, ncmx, ncmy, nthrd)
    call dmalloc(tnout, ncmx, ncmy, nthrd)
    call dmalloc(tnabsr, ncmx, ncmy, nthrd)
    call dmalloc(nfs, ncmx, ncmy)
    call dmalloc(nin, ncmx, ncmy)
    call dmalloc(nout, ncmx, ncmy)
    call dmalloc(nabsr, ncmx, ncmy)
#endif

    
!    ! accumulated arrays for GCMFD
!    call dmalloc (atrcwx,    2, ncmx,     ncmy)
!    call dmalloc (atrcwy,    2, ncmx,     ncmy)
!    call dmalloc0(aphi2x,  1,2, 0,ncmx, 0,ncmy)
!    call dmalloc0(aphi2y,  1,2, 0,ncmx, 0,ncmy)
!    call dmalloc0(aphi2trx,1,2, 0,ncmx, 0,ncmy)
!    call dmalloc0(aphi2try,1,2, 0,ncmx, 0,ncmy)
    
    ! uniform fission source
    call dmalloc(pwadj, ncmx, ncmy)
    call dmalloc(pwadjr, ncmx, ncmy)
    call dmalloc(nfsrc, ncmx, ncmy)
    call dmalloc(nfsrcth, ncmx, ncmy, nthrd)
    call dmalloc(nfsrcf, nfmx, nfmy)
    call dmalloc(nfsrcfth, nfmx, nfmy, nthrd)
    
    
    pwadj=1._8
    do j=1,ncmy
        do i=1, ncmx
            pwadjr(i,j)=1._8/pwadj(i,j)
        enddo
    enddo
    


    call dmalloc(SHNp2bc,ncmx, ncmy)
    call dmalloc(p2b,rcx, rcy, ncmx, ncmy)
    call dmalloc(powcm, ncmx, ncmy)
    call dmalloc(recpow, nfmx, nfmy)
    call dmalloc(powc, ncmx, ncmy)
    

contains

    subroutine toupper1(aa)
    !
    ! convert lowercase string to upZpercase
    integer:: INDXA, IDNXZ, lenaa, ia
    parameter (INDXA=97,IDNXZ=122)
    character aa*(*)
    !
    lenaa=len_trim(aa)
    i=1
    do while (aa(i:i).ne.' ' .and. i.le.lenaa)
     ia=ichar(aa(i:i))
     if(ia.ge.INDXA) aa(i:i)=char(ia-32)
     i=i+1
     if(i.gt.lenaa) return
    enddo
    end subroutine

    integer function nfields(aline)
    character(256)  aline
    logical nonblankd,nonblank,multidata
    integer ncol, n, multidcol, nmult
    !
    nonblankd=.false.
    multidata=.false.
    oneline=aline
    ncol=len_trim(aline)
    n=0
    do i=ncol,1,-1
     if(sc(i).eq.BLANK .or. ichar(sc(i)).eq.9) then !ichar(tab)=9
        nonblank=.false.
     else
        if(sc(i).eq.AST) then
           multidata=.true.
           multidcol=i
        endif
        nonblank=.true.
        if(sc(i).eq.BANG) then
           n=-1
           nonblankd=.true.
        endif
     endif
     if((.not.nonblankd.and.nonblank) .or. &
          (nonblankd.and..not.nonblank)) then
        n=n+1
        if(multidata.and.(nonblankd.and. .not.nonblank)) then
           read(oneline(i+1:multidcol-1),*) nmult
           n=n+(nmult-1)*2
           multidata=.false.
        endif
     endif   
     nonblankd=nonblank
    enddo
    if(mod(n,2).ne.0) then
    nfields=n/2+1
    else
    nfields=n/2
    endif
    !
    return
    end function

end subroutine
